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A8STRACT 


A  solution  procadura  Is  prasantad  for 
tha  computation  of  dynamic  stall  phenomena 
ancountarad  by  arbitrary  shapad  airfoils 
undar  arbitrary  flow  conditions.  This 
procadura  solvas  tha  unstaady,  Incomprasslbla 
Navier-Stokes  and  tha  unstaady  boundary  layar 
aquations  using  an  officiant,  zonal  approach. 
A  numbar  of  rosults  for  a  modi f lad  NACA  0012 
airfoil  axporlanclng  dynamic  stall  aro 
prasantad  and  comparad  with  aval  labia 
numarlcal  data.  Quail tatlva  comparisons  with 
flow  visualization  oxparlmants  aro  also 
prasantad.  Tha  prasant  study  also  Illustrates 
tha  effect  of  numerical  viscosity  on  tha 
accuracy  and  robustness  of  tha  solution 
procedure. 


INTRODUCTION 

The  problem  of  dynamic  stall  Is  an 
Important  area  of  research  in  tha  helicopter 
Industry,  because  of  the  largo  load 
variations,  particularly  the  pitching  moment 
variations,  that  occur  during  this 
phenomenon.  Presently,  tha  numerical  modeling 
of  this  phenomenon  Is  dona  primarily  through 
a  synthesis  of  existing  dynamic  stall 
experimental  data  .  This  approach  Is  highly 
empirical  In  nature.  It  also  relies  on  tha 
availability  of  a  large  body  of  experimental 
data  covering  a  wide  range  of  airfoil  shapes 
and  flow  conditions. 

One  of  tha  earliest  attempts  to 
numerically  simulate  tha  dynamic  stall 
phenomenon  was  made  by  Mehta  (Raf.  1).  In 
this  work  tha  Incompressible,  laminar 
Navlar-Stokas  equations  wars  solved  In  the 
vortl city-stream  function  form  using  a  finite 
difference  solution  procedure.  This  approach 
was  able  to  predict  the  major  features  of  the 
dynamic  stall  phenomenon,  Including  the 
formation  and  shedding  of  a  strong  leading 
edge  vortex.  The  laminar,  compressible 
Navlar-Stokas  equations  were  used  to  compute 
the  dynamic  stall  by  Sankar  and  Tassa  (Raf. 
2).  The  effect  of  turbulence  was  accounted 
for  In  the  works  of  Shamroth  (Ref.  3)  and 
Sankar  and  Tang  (Ref.  4).  The  above 
approaches  are  all  based  on  the  finite 
difference  solution  of  the  governing 
equations  and  require  large  amounts  of 
computer  time  and  memory  resources  to 
accurately  predict  the  dynamic  stall 
phenomenon. 
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In  the  present  work,  a  zonal  solution 
procedure  is  used  to  perform  dynamic  stall 
calculations.  This  approach  was  previously 
developed  for  stationary  airfoils  of 
arbitrary  shape  experiencing  massively 
separated  flows  (Ref.  5).  Recently  a  number 
of  improvements  have  been  made  to  this 
approach  to  reduce  the  number  of  grid  points 
and  the  computer  time  required  to  predict 
separated  flows  (Ref.  6). The  present  work 
addresses  the  extension  of  the  algorithms  and 
concepts  Implemented  In  References  5  and  6  to 
moving  grids,  and  oscillating  airfoils.  Such 
an  extension  has  resulted  In  a  robust,  useful 
solver  capable  of  generating  the  dynamic 
stall  load  hysteresis  loops  in  less  than  30 
minutes  on  a  scalar  machine  such  as  the  CDC 
CYBER  855  system. 

The  zonal  approach  used  in  this  study 
was  motivated  by  a  numbar  of  factors.  In  many 
static  and  dynamic  stall  problems,  the 
separated  flow  Is  confined  to  one  of  the 
airfoil  surfaces  (upper  or  lower)  and  the 
wake.  On  the  other  surface  the  flow  is 
attached,  and  may  be  approximated  by  an 
unsteady  boundary  layer  flow.  The  zonal 
approach  presented  here  solves  the 
computationally  costly  Navier-Stokes 
equations  only  In  the  separated  regions.  On 
the  attached  flow  side,  the  boundary  layer 
equations  are  solved.  This  reduces  the  number 
of  points  where  the  Navier-Stokes  equations 
are  solved  by  approximately  a  factor  of  two. 
Secondly,  since  in  the  separated  flow  regions 
the  length  scales  of  the  vortices  of  interest 
are  large  compared  to  the  length  scales  in 
the  boundary  layer  regions,  a  coarser  grid 
may  be  used  in  the  separated  region.  This 
also  translates  into  fewer  points  and  larger 
time  steps  in  the  separated  regions. 

A  second  feature  of  the  present  zonal 
approach  is  the  ability  to  bring  the  far 
field  boundaries  closer  to  the  airfoil.  This 
Is  achieved  by  a  closed  form  specification  of 
the  velocities  and/or  stream  function  values 
on  the  far  field  boundary  using  the 
Blot-Savart  law.  The  use  of  smaller 
computational  domains  allow  the  grid  points 
to  be  packed  close  to  the  solid  surface,  and 
also  allows  the  grid  to  stretch  smoothly  from 
the  airfoil  surface  to  the  far  field.  This  is 
an  Important  consideration  in  numerical 
solutions  if  second  order  spatial  accuracy  is 
to  be  guaranteed. 

The  zonal  approach  uses  an  integral  form  of 
the  kinematic  equations  to  determine  the 
velocity  and  stream  function  values  in  the 


separated  regions,  thus  removing  tha  naad  to 
obtain  tha  numerical  solution  of  tha 
Poisson's  aquation  for  tha  stream  function 
Iteratively.  In  tha  past,  Intagral  approachas 
such  as  tha  ona  atnployad  hara  Mara  not  as 
officiant  as  axpoctad  bacausa  of  tha  naad  to 
rapaatadly  compute  or  stora  tha  geosMtrlc 
Influence  coafflclants.  In  tha  prasant  work, 
a  Fourier  sari  as  axpanslon  of  tha  vortlclty 
flald  along  ona  of  tha  coordinate  dlractlons 
Is  usad  to  simplify  tha  aval uat ion  of  tha 
valocltlas. 

Another  Important  feature  of  tha  present 
zonal  approach  Is  tha  procedure  for  tha 
determination  of  tha  surface  vortlclty, 
generated  at  each  time  step  to  satisfy  tha  no 
slip  conditions.  The  procedure  described  hara 
ensures  that  tha  global  conservation  of 
vortlclty  Is  satisfied.  This  procedure  allows 
tha  surface  vortlclty  distribution  to  be 
determined  with  good  accuracy,  and  also 
permits  the  distribution  of  surface  pressures 
which  depend  on  tha  vortlclty  gradient  at  tha 
surface,  to  be  calculated  accurately. 

In  tha  following  sections,  tha  governing 
aquations  and  tha  numerical  procedure  are 
described.  In  a  separata  section,  the 
computer  time  and  memory  requirements  are 
given.  Finally,  a  number  of  numerical  results 
are  presented  for  the  dynamic  stall 
phenomenon  experienced  by  a  modified  NACA 
0012  airfoil,  and  compared  with  available 
numerical  data  and  flow  visualization 
studies. 

MATHEMATICAL  FORMULATION 
Governing  Equations: 

In  order  to  handle  arbitrary  airfoils 
undergoing  arbitrary  motion,  a  body-fitted, 
orthogonal  0-grld  system  Is  first 
constructed.  In  the  present  work,  the  above 
coordinate  system  was  constructed  through  an 
analytical  or  numerical  conformal  mapping  of 
the  airfoil  shape  onto  a  unit  circle, 
followed  by  a  distribution  of  nodes  on  the 
exterior  of  the  unit  circle.  This 
distribution  Is  such  that  a  sufficiently 
large  number  of  nodes  are  clustered  In  the 
vicinity  of  the  airfoil.  The 
(r,  e)  circle  plane  used  In  this  work  Is 
referred  to  as  the  c-  plane.  The  physical 
plane  (x,y)  Is  referred  to  as  the  z-  plane. 

The  governing  equations  for  the  unsteady 
Incompressible  viscous  flow  past  a  rotating 
airfoil  take  the  following  form  In  the  circle 
plane  : 

?  x  I«  v 


V  x  v  * 


H*wt  ♦  (I  X  *  W2w 


In  the  above  equations,  H  Is  the  scale  factor 
of  transformation  given  by 


H  =  (Xrya  -  x9yr)/r 


This  quantity  may  be  determined  analytically 
If  the  conformal  mapping  Is  a  Joukowski 
transformation.  It  may  be  numerically 
determined  otherwise.  The  quantity  *'  Is  the 
stream  function  In  the  rotating  frame  of 
reference,  and  Is  related  to  the  stream 
function  *  In  a  stationary  frame  of  reference 
through  the  following  relationship: 

«  ♦  ♦  Q(x2  ♦  y2)/2  (3) 


Finally,  u  Is  the  vortlclty 
distribution,  as  observed  in  an  inertial 
frame  of  reference,  and  Q  Is  the  angular 
velocity  of  the  airfoil. 


Boundary  Layer  Approximation: 

In  the  applications  considered  here,  the 
viscous  flow  region  over  the  lower  surface  is 
confined  to  a  very  thin  boundary  layer. 
Favorable  pressure  gradients  exist  over  most 
of  the  entire  lower  surface.  For  this  reason, 
In  the  present  application,  the  streamwise 
diffusion  terms  appearing  In  the  vortlclty 
transport  equation  were  neglected,  giving  the 
following  unsteady,  parabolic  partial 
differential  equation  for  u: 


H‘wt  MTxf) 


The  above  assumption  allows  the  vorticity 
field  In  the  boundary  layer  region  to  be 
determined  using  a  simple  non-iterative 
marching  scheme,  starting  from  the  leading 
edge  stagnation  point. 

Integral  Formulation  of  the  Kinematics: 

The  kinematic  relationships  for  the  stream 
function  may  be  reexpressed  as  an  integral 
expression  for  the  velocity  field  in  the 
circle  plane.  If  this  Is  done,  the  following 
Integral  equation  results: 


x  (r  -  r  ) 


lf  ’  r#l‘ 


r.dr  d6  + 
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J  r  (VxT)-yr-y  -  (VxT)x(nfl)x(r-ro) 


i?  -  y2 


In  the  above  equation,  7  Is  the  position 
vector  of  the  point  where  the  velocity  is 
computed.  Also,  r.  is  the  variable  of 
Integration.  The  term  f  contains  the 
contributions  of  the  freestream.  The  line 
integral  over  the  unit  circle  Is  zero  for 
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stationary  airfoils,  but  is  a  non-zaro 
quantity  for  oscillating  airfoils.  Th#  term 
n  appearing  in  the  above  line  integral  is  an 
outward  pointing  unit  vector. 

Integral  Formulation  for  the  Surface 
Vorticlty: 

The  surface  vorticlty  generated  at  every 
time  step  should  satisfy  the  zero  normal  and 
tangential  velocity  conditions.  The  law  of 
conservation  of  vorticlty  should  also  be 
satisfied.  In  the  present  application,  the 
Integral  relationship  for  the  velocity  at  the 
interior  points  may  be  used  directly  to 
determine  the  surface  voriticity 
distribution.  The  details  of  the  surface 
vorticlty  determination  are  given  in  the 
following  section. 


co(r)  s  j  a0  T  dro  +  co(1)/r  (8) 
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For  a  complete  list  of  the  recursive 
relationships  between  the  various  Fourier 
coefficients,  the  reader  is  referred  to 
Reference  7. 

Once,  the  coefficients  of  the  Fourier 
series  expansions  are  determined,  tne 
velocity  values  at  any  point  in  the  flow 
field  may  be  found.  In  the  present  work, 
this  discretization  was  used  at  all  the  nodes 
in  the  computational  domain  where  velocity 
values  are  needed,  including  the  nodes  within 
the  boundary  layer  region. 


NUMERICAL  FORMULATION 

Discretization  of  the  Integral  Equation  for 
Velocity: 


In  order  to  evaluate  the  velocity  field  V  at 
any  point  in  the  computational  field,  the 
following  strategy  was  used.  The  grid 
generated  using  the  mapping  procedure  My  be 
thought  of  as  a  collection  of  cells.  Inside 
every  cell  the  vorticlty  value  is  considered 
invariant  in  the  radial  direction.  In  the 
e-direction,  it  was  assumed  that  the 
vorticlty  field  My  be  approxiMted  by  the 
following  finite  Fourier  series: 


2  Vr)  r 

u(r,  8)H  *  — 2 —  +  )  [<*n(r)cosn8+@ns1n  n8] 
n«l  (6) 


The  two  components  of  the  velocity  vector 
were  also  expressed  as  the  following  Fourier 
Series  expansions  : 


«0(r)  ? 
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When  the  above  series  expansions  are 
substituted  into  Equation  (5),  the  resulting 
expressions  My  be  analytically  Integrated 
both  in  the  8-  and  the  radial  directions. 
This  yields  a  simple,  recursive  relationship 
between  the  coefficients  of  the  series  for 
the  vorticlty,  and  the  coefficients  of  the 
series  for  the  velocity  field.  For  example, 
on  all  node  points  located  on  a  circle  given 
by  r*  constant,  one  obtains. 


Determination  of  Surface  Vorticlty: 

The  surface  vorticlty  values  which  are 
needed  at  every  time  step  were  obtained  as 
follows.  The  integral  equation  for  the 
velocity  is  applied  at  the  airfoil  surface, 
given  by  r*l  in  the  circle  plane: 


r  dr  de 
0  0 


0 


I(r) 

(9) 


The  right  hand  side  of  the  above 
equations  is  the  contribution  due  to  the 
velocity  of  the  fluid  at  the  airfoil  surface 
in  the  transfoneed  plane,  and  is  entirely  due 
to  the  rotation  of  the  airfoil  about  a 
pitching  axis. 

The  Fourier  series  expansions  for  the 
vorticlty  field  are  now  substituted  into  the 
above  equation.  The  right  hand  side  of  this 
equation  is  also  expressed  as  a  finite 
Fourier  series  expansion.  All  the  terms  which 
are  multiples  of  the  ssm  sine  or  cosine  term 
are  then  grouped  together,  and  each  group  of 
terms  is  individually  set  to  zero.  If  the 
above  procedure  is  followed,  then  2N 
equations  result  for  the  2N+1  coefficients  of 
the  Fourier  series  for  the  vorticlty 
distribution  at  the  airfoil  boundary.  The 
additional  equation  required  to  determine  all 
the  2N+1  coefficients  uniquely  is  the  law  of 
total  conservation  of  vorticlty,  given  by: 


\\  Weo  *  2aA  3  0  (10) 


Where  A  is  the  area  enclosed  by  the  airfoil 
in  the  physical  plane.  This  approach  gives 
an  explicit  relationship  for  the  Fourier 
coefficients  for  the  vorticlty  distribution 
at  the  airfoil,  in  terms  of  the  vorticity 
field  away  from  the  airfoil  surface.  Note 
that  this  approach  is  consistent  with  the  way 
the  Interior  velocity  components  are 
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evaluated,  and  avoids  any  non-uni  qua 
specification  of  tha  surface  vorticity.  For  a 
set  of  explicit  relationships  between  the 
Fourier  coefficients  of  the  surface  vorticity 
and  the  interior  vorticity,  the  reader  is 
referred  to  Ref.  7. 

Numerical  treatment  of  the  Vorticity 
Transport  Equation: 

In  the  separated  flow  regions,  the 
vorticity  transport  equation  was  discretized 
as  follows: 


H2  K.  +  i  $eu>  -  i  6gp  lpw 
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In  the  above  discretization,  ,  «  etc. 
are  the  standard  central  differences'^  which 
take  into  account  the  fact  that  the  grid  is 
not  uniformly  spaced  In  the  radial  direction. 
The  quantity  ?T  is  the  backward  difference 
operator  with  respect  to  time.  The  operators 
r  and  are  four-point  upwind  difference 
operators,  patterned  after  the  QUICK  upwind 
scheme  proposed  by  Leonard  (Ref.  8). 

All  the  vorticity  values  appearing  in 
the  upwind  differences  and  the  viscous  terms 
were  kept  at  the  unknown  time  level.  Thus  the 
vorticity  values  at  all  the  nodes  are  coupled 
to  each  other  in  a  fully  implicit  manner. 
These  vorticity  values  were  iteratively 
solved  for,  using  a  successive  line  under 
relaxation  point  scheme.  The  values  of  the 
vorticity  at  the  boundary  ware  also  updated 
simultaneously  with  the  interior  vorticity 
values. 


In  tha  boundary  layer  regions,  the  same 
discretization  was  used,  except  the  diffusion 
terms  along  the  8-direction  were  suppressed. 
Since  the  flow  velocities  relative  to  the 
grid  are  always  from  the  leading  edge  to  the 
trailing  edge  on  the  lower  surface  for  the 
low  reduced  frequencies  considered  here,  the 
vorticity  transport  equation  in  the  boundary 
layer  region  may  be  solved  using  a  simple 
marching  schema.  Since  the  boundary  vorltlcty 
values  In  the  boundary  layer  region  are 
coupled  strongly  to  the  vorticity  values  in 
the  separated  region,  these  values  change 
from  one  iteration  to  another.  Thus,  It  is 
necessary  to  solve  the  boundary  layer 
equations  iteratively,  along  with  the 
vorticity  transport  equation  in  the  separated 
flow  region. 


All  the  calculations  were  carried  out  by 
starting  the  flow  from  rest  implusively,  and 
marching  in  time  until  a  steady  state  or  a 
periodic  solution  is  achieved.  In  the  case  of 
the  dynemic  stall  calculations,  a  steady 
solution  was  first  obtained  at  the  lowest 
angle  of  attack  experienced  by  the  airfoil 
during  the  dynamic  stall.  This  steady  state 


solution  was  used  as  the  initial  condition 
for  the  oscillating  airfoil  problem. 

When  advancing  the  solution  from  one 
time  step  to  the  next,  the  following 
procedure  is  followed: 

1.  The  velocity  values  at  all  the  interior 
nodes  are  computed  using  the  Fourier 
Series  expansion  approach  discussed 
earlier.  The  velocity  values  at  the  far 
field  boundary  were  also  updated.  The 
contributions  of  any  vorticity  that  had 
left  the  computational  domain  through 
convection  was  negel acted  in  this  steo. 

2.  The  vorticity  values  at  the  interior 
were  updated  (both  the  separated  and  the 
boundary  layer  region)  to  get  a  first 
estimate  of  the  vorticity  field  at  the 
new  time  level.  This  value  was 
under-relaxed  by  an  user  input  under 
relaxation  factor. 

3.  Tha  vorticity  values  at  the  solid 
surface  were  updated,  to  be  consistent 
with  the  interior  vorticity  values. 

Steps  2  and  3  are  repeated  a  number  of 
times  until  the  vorticity  values  at  all  the 
interior  nodes  and  boundary  nodes  are  fully 
converged. 

For  additional  details  of  the  solution 
procedure,  and  its  application  to  separated 
flow  problems,  the  reader  is  referred  to 
References  5  and  6. 


RESULTS  AND  DISCUSSION 

All  the  calculations  presented  here  are 
for  a  modified  NACA  0012  airfoil.  This 
airfoil  was  chosen,  because  some  well 
documented  numerical  results  for  the  laminar 
dynamic  stall  phenomenon,  and  some  water 
tunnel  flow  visualization  studies  are 
available  for  this  airfoil  (Ref.  1,9).  This 
airfoil  may  be  mapped  onto  a  unit  circle 
using  the  following  Joukowski  transformation: 


z  »  C  *  ir  ♦ 


Here  z  *  x  ♦  1  y,  and  C  =  r  exp(-  i8).  The 
constants  C,  y  etc.  determine  the  type  of  the 
resulting  airfoil  shape.  By  a  careful 
selection  of  these  coefficients,  the  airfoil 
thickness,  its  camberline  shape,  leading  and 
trailing  edge  radii  etc.  may  be  controlled. 

The  airfoil  surface  was  represented  by 
SO  nodes  in  the  computations,  clustered  near 
the  leading  and  the  trailing  edge  for  maximum 
accuracy.  In  the  radial  direction  60  nodes 
were  located.  The  stretching  in  the  radial 
direction  was  such  that  a  minimum  of  IS  nodes 
were  located  within  the  boundary  layer  on  the 
lower  surface.  The  location  of  the  first 
point  of  the  wall  in  the  circle  plane  was 
typically  between  0.006  and  0.003  units. 


rW  f"W  IT*. 


All  tne  oscillating  airfoil  calculations 
were  aone  atsout  a  mean  angle  of  attack  of  10 
degrees,  ana  an  amplitude  of  oscillation 
equal  to  10  degrees. Tne  following  cases  were 
considered: 

1.  Reynolds  NumOer  5000,  Reduced  Frequency 
based  on  chord  equal  to  0.25. 


2. 

Reynolds 

frequency. 

number 

10000, 

0.25 

reduced 

3. 

Reynolds 

frequency 

Number 

1000, 

0.25 

reduced 

4. 

Reynolds 

frequency 

number 

5000, 

0.5 

reduced 

the  airfoil  trailing  edge  as  snown  in  Figures 
22  and  2f  and  in  tne  lift  curve  (Figure  3a). 
The  pitcning  moment  recovers  nowever,  as  tne 
lift  drops,  as  seen  in  Figure  3b. 

During  tne  oownstroke,  the  flow  on  tne 
airfoil  gradually  reattacnes  from  tne  leading 
edge  to  the  trailing  edge.  The  laminar  flow 
on  the  airfoil  is  unable  to  withstand  tne 
relatively  small  adverse  pressure  gradients 
that  exist  during  this  phase  of  the  flow,  and 
a  sequence  of  small  vortices  of  positive  and 
negative  strength  are  shed  into  the  wake. 
This  accounts  for  the  large  oscillations  in 
the  lift  and  moment  forces  shown  in  Figures 
3a  and  3b.  This  is  in  contrast  to  turbulent 
flows  (Ref.  3,  4)  where  similar  variations 
are  not  observed  during  the  downstroke. 


A  steady  state  solution  at  zero  angle  of 
attack  was  supplied  as  the  initial  condition 
in  all  the  dynamic  stall  calculations. 


In  Figures  1  and  2  the  constant 
vorticity  contours  and  streamlines  in  the 
rotating  coordinate  system  are  shown  for  case 
1.  In  Figure  3,  the  integrated  airfoil  loads 
are  also  shown. These  figures  reveal  the  slow 
forward  progression  of  the  separation  point 
on  the  upper  surface  as  the  airfoil  pitches 
up.  A  large  amount  of  counterclockwise 
(positive)  vorticity  is  being  shed  into  the 
wake  during  this  phase  of  the  upstroke,  upto 
18  degrees  airfoil  incidence.  The  separated 
region  remains  confined  to  a  narrow  region  on 
the  upper  surface.  The  airfoil  lift  continues 
to  grow  during  this  phase  of  the  upstroke  as 
shown  in  Figure  3. 

Between  18  and  19  degrees  during  the 
upstroke,  the  shedding  of  the  positive 
vorticity  into  the  wake  stops,  and  a  region 
of  counterclockewise  vorticity  begins  to  grow 
on  the  upper  surface  as  seen  In  Figures  lb 
and  2C.  The  airfoil  lift  coefficient 
experiences  a  momentary  drop  as  shown  In 
Figure  3a. 

Subsequently,  a  leading  edge  vortex  of 
clockwise  sense  begins  to  grow  near  the 
quarter  chord  as  shown  In  Figure  lb.  The 
strength  of  the  vortex  increases  with  the 
airfoil  incidence.  The  lift  coefficient 
begins  to  grow,  but  the  moment  coefficient 
remains  relatively  unaffected.  This  Is 
because  the  pressure  suction  peaks  associated 
with  the  vortex  are  too  close  to  the  moment 
axis  (quarter  chord)  to  Influence  the  moment 
coefficient. 

As  the  leading  edge  vortex  grows  In 
strength  It  begins  to  drift  downstream,  at  a 
speed  equal  to  approxlmatley  half  the 
freestream  velocity.  The  lift  coefficient 
continues  to  grow  as  shown  In  portion  BC  of 
the  lift  curve  In  Figure  3a.  The  moment 
coefficient  begins  to  Increase  In  magnitude 
because  the  suction  peaks  associated  with  the 
vortex  are  now  sufficiently  far  away  from  the 
pitching  axis,  and  eventually  moment  stall 
occurs,  as  shown  in  portion  BC  In  Figure  3b. 
The  lift  stall  occurs  when  the  vortex  reaches 


The  drag  variations  during  the  dynamic 
stall  process  are  shown  in  Figure  3c.  Except 
at  small  angles  of  attack,  the  primary 
contribution  to  the  drag  is  from  the  pressure 
drag. 

Some  numerical  experiments  were  done  to 
determine  the  effect  of  numerical  viscosity 
due  to  the  upwind  differencing  of  the 
convection  terms,  on  the  solution.  In  Figure 
4,  the  load  hysteresis  for  case  1,  computed 
with  a  first  order  accurate  upwind  difference 
scheme  is  shown.  The  large  numerical 
viscosity  associated  with  the  first  order 
scheme  smears  out  a  number  of  features 
observed  using  the  third  order  accurate 
upwind  scheme,  and  only  a  qualitative 
resemblance  between  Figures  3  and  4  is  found. 

Case  2,  is  similar  to  case  1  except  the 
flow  Reynolds  number  is  increased  from  5000 
to  10000.  In  Figure  5  the  lift  and  moment 
hysteresis  loops  are  shown  for  case  2. 

Case  3  Is  similar  to  cases  1  and  2 
except  the  Reynolds  number  was  1000.  At  this 
lower  Reynolds  number  no  clearly 
distinguishable  leading  edge  vortex  structure 
was  observed,  as  seen  In  the  streamline 
contours  (Figure  6).  The  associated  lift  and 
moment  loops  are  plotted  In  Figure  7. 

This  type  of  dynamic  stall  is  often 
classified  as  the  trailing  edge  stall. 

The  last  set  of  dynamic  stall 
calculations  presented  here  are  for  the  same 
flow  conditions  as  case  1,  except  at  the 
slightly  higher  reduced  frequency  of  0.5. 
This  case  was  chosen  because  of  the 
availability  of  numerical  data  (Ref.  1)  and 
water  tunnel  flow  visualization  pictures.  The 
physical  phenomena  at  this  reduced  frequency 
are  very  similar  to  case  1.  The  actual  drop 
in  lift  and  moment  coefficients  are  however 
smaller  than  at  the  lower  reduced  frequency. 
In  Figure  8,  the  lift,  drag  and  moment 
hysteresis  loops  are  plotted  and  compared 
with  the  results  of  Mehta.  Good  quantitative 
agreement  is  observed  during  the  upstroke, 
but  some  discrepancies  between  the  two  set  of 
data  Is  found  during  the  downstroke.  In  view 
of  the  differences  In  the  solution  procedures 
used  in  Reference  1  and  in  this  work,  these 
differences  are  not  unreasonable.  In  Figure 
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9,  the  streamlines  in  an  inertial  frame  are 
plotted  at  selected  time  levels  and  compared 
with  the  water  tunnel  observations  of  Werle' 
presented  in  Reference  1.  Excellent  agreement 
between  these  two  sets  of  data  is  found. 

In  all  the  cases  computed  above, 
computer  time  per  case  ranged  between  15  and 
30  minutes  on  a  COC  CYBER  855,  depending  on 
the  reduced  frequency  and  Reynolds  number. 
This  code  is  now  being  vectorized  for  a  COC 
CYBER  205  computer  system,  and  preliminary 
calculations  show  that  the  computer  time  per 
cycle  can  be  reduced  to  less  than  2  minutes 
on  this  system. 


CONCLUDING  REMARKS 

An  efficient  solution  procedure  has  been 
developed  for  the  prediction  of  the  dynamic 
stall  characteristics  of  arbitrary  airfoils 
in  laminar,  incompressible  flow.  This 
solution  procedure  is  an  order  of  magnitude 
more  efficient  than  existing  finite 
difference  procedures  for  the  same  problem.  A 
two  layer  eddy  viscosity  model  is  being 
implemented  into  this  solver  and  will  be  used 
to  study  turbulent  dynamic  stall. 

The  computations  reveal  characteristics 
similar  to  turbulent  leading  edge  stall  at 
Reynolds  numbers  5000  and  10000,  while  at  the 
lower  Reynolds  number  of  1000,  a  trailing 
edge  stall  phenomenon  is  observed.  The  effect 
of  higher  reduced  frequencies  Is  to  reduce 
the  lift  and  moment  drops  when  stall  occurs. 
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(a)  Angle  of  Attack  =  16.90' 


(b)  Angle  of  Attack  =  18.34> 


(d)  Angle  of  Attack  =  19.53. 


Figure  1.  Vorticity  Contours  around  an  OscMating 
NACA  0012  Airfoil.  Reynolds  Numoer  5000,  0.25 
Reduced  F-equency. 


F1gur«  5.  Lift  and  Moment  Hysteresis  Loops  at 
Reynolds  Number  10000,  0.25  Reduced  Frequency. 
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Figure  7.  Lift  and  Moment  Hysteresis  loops  for  a 
NACA  0012  Airfoil.  Reynolds  Number  1000,  0.25 
Reduced  Frequency. 
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Figure  9.  Comparison  of  Streamline  in  the  Inertial 
Frame  with  Water  Tunnel  Observations.  Reynolds 
Number  5000,  0.5  Reduced  Frequency. 
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Figure  8.  Comparison  of  Lift  and  Moment  Hysteresis 
Loops  Computed  Using  the  Present  Scheme  with  the 
Results  of  Mehta.  Reynolds  Number  5000,  0.5 
Reduced  Frequency. 
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